rm(list = ls())
require(foreign)
require(stats4)
require(mfx)
require(stargazer)

#import data
symmetrydata=read.dta("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Expt1.dta")
symmetryframe=data.frame(id=symmetrydata$var1,user_id=symmetrydata$var2,qn=symmetrydata$var3,qid=symmetrydata$var4,order=symmetrydata$var5,bid=symmetrydata$var6,chosen_act=symmetrydata$var9, true_state=symmetrydata$var10)

#reward value
r=10 

#data frame for letters
symframeletters=subset(symmetryframe,user_id<1700)
symframeletters$correct=(symframeletters$true_state<=13 & symframeletters$chosen_act==1)+(symframeletters$true_state>13 & symframeletters$chosen_act==2)
unique(symframeletters$qid)

#data frame for balls
symframeballs=subset(symmetryframe,user_id>1700)
symframeballs$true_state=symframeballs$true_state+6*(symframeballs$qid==3)+4*(symframeballs$qid==4)+2*(symframeballs$qid==5)
symframeballs$correct=(symframeballs$true_state<=10 & symframeballs$chosen_act==1)+(symframeballs$true_state>10 & symframeballs$chosen_act==2)
ballsqid=unique(symframeballs$qid)

#Variables to be stored
acculin=matrix(rep(0,80),nrow=4,ncol=20)
accupow=matrix(rep(0,80),nrow=4,ncol=20)
accugen=matrix(rep(0,80),nrow=4,ncol=20)
accuexp=matrix(rep(0,80),nrow=4,ncol=20)
accunoneigh=matrix(rep(0,80),nrow=4,ncol=20)
dataaccu=matrix(rep(0,80),nrow=4,ncol=20)


#Qids for each treatment type
lettersqid=c(3,5,7,8)
ballsqid=c(3,4,5,6)

#number of states used by different treatments
Nvec=c(8,12,16,20)

#States for different treatments
states=list()
states[[1]]=c(7:14)
states[[2]]=c(5:16)
states[[3]]=c(3:18)
states[[4]]=c(1:20)

correct=list()
incorrect=list()

symframeballs$distance=abs(symframeballs$true_state-10.5)

#{*} Make correct vecs (needs changing between treatment types)
for(i in 1:4){
correct[[i]]=vector('numeric')
incorrect[[i]]=vector('numeric')
for(j in 1:Nvec[i]){
correct[[i]][j]=sum(symframeballs$correct*(symframeballs$qid==ballsqid[i])*(symframeballs$true_state==states[[i]][j]))
incorrect[[i]][j]=sum((1-symframeballs$correct)*(symframeballs$qid==ballsqid[i])*(symframeballs$true_state==states[[i]][j]))
}
}


multinom=function(vec){
return(sum(log(1:sum(vec)))-sum(log(c(1:vec[1],1:vec[2]))))
}
multinomlist=list()
multinomvec=vector("numeric")
for(i in 1:4){
multinomlist[[i]]=vector("numeric")
for(j in 1:Nvec[i]){
multinomlist[[i]][j]=multinom(c(correct[[i]][j],incorrect[[i]][j]))
}
multinomvec[i]=sum((multinomlist[[i]]))
}

multinomconst=sum(multinomvec)
######################
#Linear Mutual Information Costs

accuracylin=function(kappaall,kappaloc,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))
return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

likelihoodlin=function(kappaall,kappaloc){
accuvec1=accuracylin(kappaall,kappaloc,1)
accuvec2=accuracylin(kappaall,kappaloc,2)
accuvec3=accuracylin(kappaall,kappaloc,3)
accuvec4=accuracylin(kappaall,kappaloc,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 
modellin=mle(likelihoodlin,start=list(kappaall=5.1,kappaloc=5.1), method="L-BFGS-B",lower=c(0.1,0.1),upper=c(100,1000))

kappaalllin=coef(modellin)[1]
kappaloclin=coef(modellin)[2]
for(expi in 1:4){
N=Nvec[expi]
accuracy=accuracylin(coef(modellin)[1],coef(modellin)[2],expi)
acculin[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}

#################################################

#power costs

accuracypow=function(kappaall,kappaloc,rho,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
overallcost=sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))-N*1/N*log(1/N)
neighborcost=vector('numeric')
for(j in 1:(N/2-1)){
neighborcost[j]=1/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))-2*0.5*log(0.5)
}
neighborcost[N/2]=(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))-2*0.5*log(0.5)
return(kappaall*overallcost^rho+kappaloc*sum(neighborcost^rho))
}

#Optimization target
negpayoff=function(accu){

return(-1*(r*sum(accu/(N/2))-cost(accu)))

}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2)-c(1:(N/2))/200,negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

likelihoodpow=function(kappaall,kappaloc,rho){
accuvec1=accuracypow(kappaall,kappaloc,rho,1)
accuvec2=accuracypow(kappaall,kappaloc,rho,2)
accuvec3=accuracypow(kappaall,kappaloc,rho,3)
accuvec4=accuracypow(kappaall,kappaloc,rho,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

modelpow=mle(likelihoodpow,start=list(kappaall=5,kappaloc=5.5, rho=1.1), method="L-BFGS-B",lower=c(0.1,0.1,0.1),upper=c(1000,1000,100))
summary(modelpow)

kappaallpow=coef(modelpow)[1]
kappalocpow=coef(modelpow)[2]
rhopow=coef(modelpow)[3]
for(expi in 1:4){
N=Nvec[expi]
accuracy=accuracypow(coef(modelpow)[1],coef(modelpow)[2],coef(modelpow)[3],expi)
accupow[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}


#################################################
#Gen entropy approach
accuracygen=function(kappaall,kappaloc,q,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Define constraint matrix
lastmat=cbind(diag(N/2-1), matrix(0, N/2-1,1))-0.99999*cbind( matrix(0, N/2-1,1),diag(N/2-1))
conmat=rbind(diag(N/2),-diag(N/2),lastmat)
convec=c(rep(0.5,N/2),rep(-0.99,N/2),rep(0,N/2-1))

#Information cost function
cost=function(accu){
neighborcost=vector('numeric')

if(q==1){

overallcost=kappaall*sum(accu/(N/2)*log(accu/(N/2))+(1-accu)/(N/2)*log((1-accu)/(N/2)))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc/2*(accu[j]*log(accu[j]/(accu[j]+accu[j+1]))+(1-accu[j])*log((1-accu[j])/(2-accu[j]-accu[j+1]))+accu[j+1]*log(accu[j+1]/(accu[j]+accu[j+1]))+(1-accu[j+1])*log((1-accu[j+1])/(2-accu[j]-accu[j+1])))
}
neighborcost[N/2]=kappaloc*(accu[N/2]*log(accu[N/2])+(1-accu[N/2])*log(1-accu[N/2]))

}else{

if(q==2){
overallcost=kappaall*(-1/N*sum(log(accu/(N/2))+log((1-accu)/(N/2)))-2*log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(-1/2*((accu[j]+accu[j+1])/2*(log(accu[j]/(accu[j]+accu[j+1]))+log(accu[j+1]/(accu[j]+accu[j+1])))+(2-accu[j]-accu[j+1])/2*(log((1-accu[j])/(2-accu[j]-accu[j+1]))+log((1-accu[j+1])/(2-accu[j]-accu[j+1]) ) )) -2*log(2))
}
neighborcost[N/2]=kappaloc*(-1/2*(log(accu[N/2])+log(1-accu[N/2]))-2*log(2))


}else{
overallcost=kappaall*(N^(1-q)/((q-2)*(q-1))*(sum((accu/(N/2))^(2-q)+((1-accu)/(N/2))^(2-q)))-1/((q-2)*(q-1))-log(N))
for(j in 1:(N/2-1)){
neighborcost[j]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[j]+accu[j+1])/2*((accu[j]/(accu[j]+accu[j+1]))^(2-q)+(accu[j+1]/(accu[j]+accu[j+1]))^(2-q))+(2-accu[j]-accu[j+1])/2*(((1-accu[j])/(2-accu[j]-accu[j+1]))^(2-q)+((1-accu[j+1])/(2-accu[j]-accu[j+1]))^(2-q)))-1/((q-2)*(q-1))-log(2))

}
neighborcost[N/2]=kappaloc*(2^(1-q)/((q-2)*(q-1))*((accu[N/2])^(2-q)+(1-accu[N/2])^(2-q))-1/((q-2)*(q-1))-log(2))

}

}

return(overallcost+sum(neighborcost))
}

#Optimization target
negpayoff=function(accu){
return(-1*(r*sum(accu/(N/2))-cost(accu)))
}

#Find optimal accuracies
optimout=constrOptim(rep(0.6,N/2),negpayoff, grad=NULL,ui=conmat,ci=convec)
accuvec=optimout$par
return(accuvec)
}

likelihoodgen=function(kappaall,kappaloc,q){
accuvec1=accuracygen(kappaall,kappaloc,q,1)
accuvec2=accuracygen(kappaall,kappaloc,q,2)
accuvec3=accuracygen(kappaall,kappaloc,q,3)
accuvec4=accuracygen(kappaall,kappaloc,q,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}

modelgen=mle(likelihoodgen,start=list(kappaall=5,kappaloc=5, q=1), method="L-BFGS-B",lower=c(0.01,0.01,0.01),upper=c(100,100,100))

kappaallgen=coef(modelgen)[1]
kappalocgen=coef(modelgen)[2]
qgen=coef(modelgen)[3]

for(expi in 1:4){
N=Nvec[expi]
accuracy=accuracygen(coef(modelgen)[1],coef(modelgen)[2],coef(modelgen)[3],expi)
accugen[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}

######################
#Linear Mutual Information Costs no neighborhoods

accuracynoneigh=function(kappaall,treatment){
#number of states in the treatment
N=Nvec[treatment]

#Information cost function
cost=function(accu){
overallcost=kappaall*(accu*log(accu/(N/2))+(1-accu)*log((1-accu)/(N/2)))
return(overallcost)
}

#Optimization target
negpayoff=function(accuscalar){
return(-1*(r*sum(rep(accuscalar,N/2)/(N/2))-cost(accuscalar)))
}

#Find optimal accuracies
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec=rep(optimout$minimum,N/2)
return(accuvec)
}

likelihoodnoneigh=function(kappaall){
accuvec1=accuracynoneigh(kappaall,1)
accuvec2=accuracynoneigh(kappaall,2)
accuvec3=accuracynoneigh(kappaall,3)
accuvec4=accuracynoneigh(kappaall,4)

#reflect accuracy vector
accuvec1=c(accuvec1,rev(accuvec1))
accuvec2=c(accuvec2,rev(accuvec2))
accuvec3=c(accuvec3,rev(accuvec3))
accuvec4=c(accuvec4,rev(accuvec4))

#find likelihood of observed data
like=multinomconst+sum(correct[[1]]* log(accuvec1)+incorrect[[1]]*log(1-accuvec1))+sum(correct[[2]]* log(accuvec2)+incorrect[[2]]*log(1-accuvec2))+sum(correct[[3]]* log(accuvec3)+incorrect[[3]]*log(1-accuvec3))+sum(correct[[4]]* log(accuvec4)+incorrect[[4]]*log(1-accuvec4))
return(-like)
}
 

modelnoneigh=mle(likelihoodnoneigh,start=list(kappaall=30), method="L-BFGS-B",lower=1,upper=100)

kappaallnoneigh=coef(modelnoneigh)[1]

for(expi in 1:4){
N=Nvec[expi]
accuracy=accuracynoneigh(coef(modelnoneigh)[1],expi)
accunoneigh[expi,]=c(rep(0,(20-N)/2),accuracy,rev(accuracy),rep(0,(20-N)/2))
}



###################################################################
#Create the actual data with bootstrap standard errors
#Store upper and lower bounds of 95% CI
dataaccuupper=matrix(rep(0,80),nrow=4,ncol=20)
dataacculower=matrix(rep(0,80),nrow=4,ncol=20)

#number of bootstrap iterations
rep=1000

for(expi in 1:4){
#Setup for speicifc treatment
playframe=subset(symframeballs,qid=ballsqid[expi])
N=Nvec[expi]

#Define Players
players=unique(playframe$user_id)
nplay=length(players)

#Set up value storage
acculist=list()
for(i in 1:N){
acculist[[i]]=vector("numeric")
}

#Draws
drawvec=sample(c(1:nplay),nplay*rep,replace=TRUE)

for(i in 1:rep){
#define bootframe
bootframe=data.frame(id=vector("numeric"), user_id=vector("numeric"), qn=vector("numeric"), order=vector("numeric"), bid=vector("numeric"),chosen_act=vector("numeric"), true_state=vector("numeric"), correct=vector("numeric"))

#fill bootframe
for(j in 1:nplay){
bootframe=rbind(bootframe, subset(playframe, user_id==players[drawvec[(i-1)*nplay+j]] ) )
}

#caluculate accuracies
for(j in 1:N){
bootsubframe=subset(bootframe,true_state==states[[expi]][j])
acculist[[j]][i]=sum(bootsubframe$correct)/length(bootsubframe$correct)
}

}

accuracyupper=vector("numeric")
accuracylower=vector("numeric")

#sort vectors
for(i in 1:N){
acculist[[i]]=sort(acculist[[i]])
accuracyupper[i]=acculist[[i]][975]
accuracylower[i]=acculist[[i]][25]
}

#calculate
accuracy=correct[[expi]]/(correct[[expi]]+incorrect[[expi]])

#fill in accuvec
dataaccu[expi,]=c(rep(0,(20-N)/2),accuracy,rep(0,(20-N)/2))
dataaccuupper[expi,]=c(rep(0,(20-N)/2),accuracyupper,rep(0,(20-N)/2))
dataacculower[expi,]=c(rep(0,(20-N)/2),accuracylower,rep(0,(20-N)/2))
}

####################################################################################################################
barplot(dataaccu, beside=TRUE, main="Accuracy By State", xlab="True State ", ylab="Accuracy ", ylim=c(0.5,1),names.arg=c(1:20), col=c("purple","orange", "red", "cyan"),xpd = FALSE)
legend("topleft",legend=c("T 1","T 2","T 3","T 4"),lwd=c(10,10,10,10),col=c("purple","orange", "red", "blue"))
#########################################################################################################################

#treatment 1
dispmat=rbind(acculin[1,7:14],accupow[1,7:14],accuexp[1,7:14],accugen[1,7:14],accunoneigh[1,7:14], dataaccu[1,7:14])
#*****************************************
barplot(dispmat, beside=TRUE, main="Accuracy By State", xlab="True State ", ylab="Accuracy ", ylim=c(0.5,1),names.arg=c(7:14), col=c("purple","orange", "red", "blue", "cyan","black"),xpd = FALSE)
legend("top",legend=c("Mutual Info","Power","Exponential","General Entropy","No Neighbor", "Data"),lwd=c(10,10,10,10,10,10),col=c("purple","orange", "red", "blue", "cyan","black"))

#treatment 2
dispmat=rbind(acculin[2,5:16],accupow[2,5:16],accuexp[2,5:16],accugen[2,5:16],accunoneigh[2,5:16], dataaccu[2,5:16])
#*****************************************
barplot(dispmat, beside=TRUE, main="Accuracy By State", xlab="True State ", ylab="Accuracy ", ylim=c(0.5,1),names.arg=c(5:16), col=c("purple","orange", "red", "blue", "cyan","black"),xpd = FALSE)
legend("top",legend=c("Mutual Info","Power","Exponential","General Entropy","No Neighbor", "Data"),lwd=c(10,10,10,10,10,10),col=c("purple","orange", "red", "blue", "cyan","black"))

#treatment 3
dispmat=rbind(acculin[3,3:18],accupow[3,3:18],accuexp[3,3:18],accugen[3,3:18],accunoneigh[3,3:18], dataaccu[3,3:18])
#*****************************************
barplot(dispmat, beside=TRUE, main="Accuracy By State", xlab="True State ", ylab="Accuracy ", ylim=c(0.5,1),names.arg=c(3:18), col=c("purple","orange", "red", "blue", "cyan","black"),xpd = FALSE)
legend("top",legend=c("Mutual Info","Power","Exponential","General Entropy","No Neighbor", "Data"),lwd=c(10,10,10,10,10,10),col=c("purple","orange", "red", "blue", "cyan","black"))

#treatment 4
dispmat=rbind(acculin[4,1:20],accupow[4,1:20],accuexp[4,1:20],accugen[4,1:20],accunoneigh[4,1:20], dataaccu[4,1:20])
#*****************************************
barplot(dispmat, beside=TRUE, main="Accuracy By State", xlab="True State ", ylab="Accuracy ", ylim=c(0.5,1),names.arg=c(1:20), col=c("purple","orange", "red", "blue", "cyan","black"),xpd = FALSE)
legend("top",legend=c("Mutual Info","Power","Exponential","General Entropy","No Neighbor", "Data"),lwd=c(10,10,10,10,10,10),col=c("purple","orange", "red", "blue", "cyan","black"))

#Linear Mutual Information Outputs
#Power Mutual Information Outputs
#General Entropy Outputs
#No Neighbor

#some output tables
dp11frame=data.frame(state=states[[1]],data=dataaccu[1,states[[1]]],datahigh=dataaccuupper[1,states[[1]]],datalow=dataacculower[1,states[[1]]],mutual_info=acculin[1,states[[1]]],power_info=accupow[1,states[[1]]],gen_ent=accugen[1,states[[1]]],no_neigh=accunoneigh[1,states[[1]]])
stargazer(dp11frame,summary=FALSE)
write.table(dp11frame, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/ballsfits.csv", sep = ",", col.names = T, append = F)

dp12frame=data.frame(state=states[[2]],data=dataaccu[2,states[[2]]],datahigh=dataaccuupper[2,states[[2]]],datalow=dataacculower[2,states[[2]]],mutual_info=acculin[2,states[[2]]],power_info=accupow[2,states[[2]]],gen_ent=accugen[2,states[[2]]],no_neigh=accunoneigh[2,states[[2]]])
stargazer(dp12frame,summary=FALSE)
write.table(dp12frame, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/ballsfits.csv", sep = ",", col.names = T, append = T)

dp13frame=data.frame(state=states[[3]],data=dataaccu[3,states[[3]]],datahigh=dataaccuupper[3,states[[3]]],datalow=dataacculower[3,states[[3]]],mutual_info=acculin[3,states[[3]]],power_info=accupow[3,states[[3]]],gen_ent=accugen[3,states[[3]]],no_neigh=accunoneigh[3,states[[3]]])
stargazer(dp13frame,summary=FALSE)
write.table(dp13frame, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/ballsfits.csv", sep = ",", col.names = T, append = T)

dp14frame=data.frame(state=states[[4]],data=dataaccu[4,states[[4]]],datahigh=dataaccuupper[4,states[[4]]],datalow=dataacculower[4,states[[4]]],mutual_info=acculin[4,states[[4]]],power_info=accupow[4,states[[4]]],gen_ent=accugen[4,states[[4]]],no_neigh=accunoneigh[4,states[[4]]])
stargazer(dp14frame,summary=FALSE)
write.table(dp14frame, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/ballsfits.csv", sep = ",", col.names = T, append = T)


